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An algorithm to find a graded Projected Entangled-Pair State representation of the ground state wave func- 
tions is developed for translationally invariant strongly correlated electronic systems on infinite-size lattices in 
two spatial dimensions. It is tested for the two-dimensional t — J model at and away from half filling, with 
truncation dimensions up to 6. We are able to locate a line of phase separation, which qualitatively agrees with 
the results based on the high-temperature expansions. We find that the model exhibits an extended s-wave su- 
perconductivity for J = OAt at quarter filling. However, we emphasize that the currently accessible truncation 
dimensions are not large enough, so it is necessary to incorporate the symmetry of the system into the algorithm, 
in order to achieve results with higher precision. 
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PACS numbers: 02.70. -c,71.10.Fd,71.10.Pm 

The investigation of models of strongly correlated electrons 
in two spatial dimensions remains to be a major challeng- 
ing issue in condensed matter physics. Actually, no well- 
controlled analytical techniques are available to study even the 
ground state properties, which has led numerous theorists to 
appeal to numerical simulations. Up to now, some powerful 
numerical approaches to classically simulate quantum many- 
body lattice systems have been proposed, such as Quantum 
Monte Carlo (QMC) [1] and the Density Matrix Renormaliza- 
tion Group (DMRG) [2]. However, the QMC suffers from a 
notorious sign problem for both strongly correlated electronic 
systems and frustrated spin systems, whereas the DMRG is 
not so efficient for quantum lattice many-body systems in two 
spatial dimensions. 

Recently, significant advances have been made in the con- 
text of classical simulations of quantum lattice many-body 
systems in terms of the so-called Tensor Network (TN) al- 
gorithms Jl fl S H 0, S EJSl- These include the Ma- 
trix Product States (MPS) 111 111 for quantum lattice systems 
in one spatial dimension, the Projected Entangled-Pair States 
(PEPS) [6] for quantum lattice systems in two and higher spa- 
tial dimensions, and the Multi-scale Entanglement Renormal- 
ization Ansatz (MERA) [10] for quantum lattice systems in 
any spatial dimensions. One of the advantages of the TN algo- 
rithms is that, in contrast to the QMC, they do not suffer from 
any sign problem, although a graded version of the TN algo- 
rithms is necessary to take into account all signs arising from 
the anti-commutivity of fermionic operators at different lat- 
tice sites. Therefore, it is highly desirable to develop efficient 
graded TN algorithms that enable us to classically simulate 
quantum electronic lattice systems in two spatial dimensions. 
Remarkably, algorithms to tackle signs arising from the anti- 
commutivity of fermionic operators at different lattice sites 
for strongly correlated electronic systems have recently been 
proposed in the context of the MERA representations Il2ll . 

In this paper, we develop a numerical algorithm to find a 
graded Projected Entangled-Pair State (gPEPS) representation 
of the ground state wave functions for translationally invariant 



strongly correlated electronic systems on infinite-size lattices 
in two spatial dimensions. In our opinion, the gPEPS is a nat- 
ural extension of the PEPS to tackle quantum electronic lat- 
tice systems, in which a parity is attached to each of the basis 
vectors of both auxiliary and physical spaces that are super 
spaces in mathematics. The algorithm is tested for the two- 
dimensional t - J model at and away from half filling, with 
truncation dimensions up to 6. We are able to locate a line 
of phase separation (PS), which qualitatively agrees with the 
results based on the high-temperature expansions IU3I1 . We 
find that the model exhibits an extended .s-wave superconduc- 
tivity for J = 0.4f at quarter filling. However, we empha- 
size that the currently accessible truncation dimensions are not 
large enough, so it is necessary to incorporate the symmetry 
of the system into the algorithm, in order to achieve results 
with higher precision. 

Graded PEPS representations. Consider a translationally 
invariant quantum electronic system on an infinite-size square 
lattice in two spatial dimensions. Suppose it consists of the 
nearest-neighbor interactions, characterized by a Hamiltonian 
H = 2<i7> h<u>- Our purpose is to find the ground state wave 
function via an imaginary time evolution, with a randomly 
chosen state as an initial state |t/<o): 



Wr) = 



exp(-//T)|l/r ) 

|exp(-ffr)hfo>l 



(1) 



when r — > oo, as long as the initial state is not orthogonal to 
the genuine ground state. 

In order to carry out the imaginary time evolution effi- 
ciently, we need to represent the system's ground state wave 
functions in terms of graded PEPS states for translationally 
invariant strongly correlated electron systems on infinite-size 
square lattices in two spatial dimensions. At each site, there 
is a local d-dimensional Hilbert super space V whose basis 
vectors are \s) {s = 1,2, • ■ • , d), with the parity [s] being 
for even vectors and 1 for odd vectors. In t he g raded version 
of the valence bond state (gVBS) picture 01411 . one may as- 
sociate four D-dimensional auxiliary super spaces V/, V r , V u , 



2 



and Vd to the physical Hilbert super space V. Suppose |Z), |r), 
\u), and \d) are bases of the auxiliary super spaces V/, V r , V u , 
and Vd, with their corresponding parities [/], [r], [u] and [d], 
respectively. Following Ref. jfl], we define a gVBS state 



(2) 



where \<f>) is a maximally entangled state \<f>) = £„ =1 
the tensor product ® is over all possible bonds on the square 
lattice, and Pis a projection operator P: V/® V,® V„®V</ — » V, 
defined as 



(3) 



l,y,u,d-l s=l 



For convenience, we assume that Wf. = if [s] + [I] + [r] + 
[u] + [d] + mod 2. Substituting Eq. © into Eq.©, and tak- 
ing into account the signs arising from the grading structure, 
under the convention that physical states | s) on a square lattice 
are arranged by first ordering from left to right along horizon- 
tal bonds and then from up to down along vertical bonds, we 
may map a gVBS to a gPEPS described by a seven-index ten- 
sor W; 



lrud\V r' ' 



w, 



lrml:t 



v = (-D 



r]([«]+[d]) 



(-l)'' r W S lrud 8t> +r >+[ u ]+[d] mod2, 0- (4) 



Here, /' and r' (l',r' = 0,1) are indices labeling two ex- 
tra horizontal grading bonds attached to each lattice site (see 
Fig- Ei))- The gPEPS for this convention is visualized in 
Fig.[TJn)- However, there exists another equivalent represen- 
tation 

Kud;,> r > = (-l) W<[ " ]+[rf]) (-1)"'' Wl ud W +M+M m od2, 0- (5) 

We emphasize that, as we shall see later on, this convention 
is only useful to absorb a two-site gate acting on a horizontal 
bond during the imaginary time evolution. In order to absorb 
a two-site gate acting on a horizontal bond during the imagi- 
nary time evolution, we need another convention that physical 
states \s) on a square lattice are arranged by first ordering from 
up to down along horizontal bonds and then from left to right 
along vertical bonds, which yields other two equivalent repre- 
sentations: 



M(M+M) 



(-1)"' W S udlr S u > + d'+[l]+[r] mod2, 0) 



and 



Kd, r yd< = (-D 



KKH+M) 



(-1)'" Kd,r 6u 



udlr °u'+d'+[t\+[r] mod2 



(6) 



o- (7) 



Here, u' and d' (u',d' = 0, 1) are indices labeling two ex- 
tra vertical grading bonds attached to each lattice site, as 
shown in Fig. [Tfiv)- Note that W* dl is related to Wf d via 



W udh = (-l) W]+[rM " ]+[d]) W* rud . The gPEPS for this conven- 
tion is visualized in Fig.[Tfv). 



Note that Eq. (0]l has been introduced in Ref. 11511 in the 
context of a fermionic PEPS (fPEPS) representation. How- 
ever, an essential difference between an fPEPS and a gPEPS 



lies in the fact that the latter may be used to absorb a two- 
site gate during the imaginary time evolution which acts on a 
horizontal bond or vertical bond (see below). In addition, it is 
convenient to use super spaces that naturally describe physical 
Hilbert spaces in the two-dimensional t — J model. 

The algorithm. As usual, the imaginary time evolution op- 
erator exp(-//r) in Eq. (Q} is implemented by dividing t into 
M small time slices St: t = MSt . For each small time 
slice St, it is represented by exp(-HST). In fact, for our pur- 
pose, we shall choose a plaquette as a unit cell, with its ver- 
tices labeled as W, X, Y and Z (see Fig. [TJiii) and Fig. QJvi))- 
Then, as follows from the Suzuki-Trotter decomposition 111 611 . 
exp(-H6r) is a product of eight different kinds of two-site 
gates U a (a = WX, XW YZ, ZY WY, YW, XZ, ZX) cor- 
responding to eight different kinds of bonds, with the two-site 
gate U a defined by 



U a = exp(— haSr), St <k 1. 



(8) 



Thus, we have reduced the problem to implement the imag- 
inary time evolution to how to update the gPEPS tensors W, 
X, Y, and Z under the action of a two-site gate U„ acting on 
eight different types of bonds. An efficient (but not optimal) 
way to do this is to adapt the strategy used in the iMPS al- 
gorithm [7]. Therefore, we attach a diagonal singular value 
matrix A a to each type of bonds, with tensors Tw, Yx, Fy, and 
Tz defined via removing a square root of the singular value 
matrix from each of all four bonds surrounding W, X, Y, and 
Z, respectively. As such, the algorithm consists of two parts: 
first, absorb the action of a two-site gate U a on a gPEPS to 
update the gPEPS tensors; second, read out the expectation 
value of a physical observable in a given gPEPS. 

(i) Updating of the gPEPS tensors. Our choice of the unit 
cell in the gPEPS representation assumes that it is transla- 
tionally invariant under two-site shifts, which implies that one 
only needs to address two consecutive sites linked by a certain 
kind of bonds; once this is done, we simultaneously update all 
the tensors on the sites linked by the same kind of bonds. The 
updating procedure for a two-site gate acting on a WX bond is 
visualized in Fig. [2] which consists of a few steps: (i) the two- 
site gate U a is applied onto the gPEPS. (ii) A single tensor © 
is formed by contracting the tensors Tw, Fx, A xw , A zx , A xz , A wy , 
A yw , and the gate U a - (iii) Reshape the tensor into a matrix 
M. (iv) A singular value decomposition (SVD) is performed 
for the matrix M, followed by a truncation, with only the D 
largest singular values retained in the updated singular matrix 
A wx . (v) Reshape the matrices U and V into the tensors U and 
V. (vi) Recover the diagonal matrix A XK , A zx , A xz , A wy , A yw , and 
update the tensors Fw and Fx to be F^ and F x . 

(ii) Measuring a physical observable. Once a gPEPS is gen- 
erated as a ground state wave function of a translationally in- 
variant quantum electronic system on an infinite-size square 
lattice, we need to compute the expectation value of a phys- 
ical observable. For this purpose, the basic building blocks 
are double tensors w, x, y and z formed from contracting the 
physical indices for the gPEPS tensors W, X, Y, and Z and 
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FIG. 1: (color online) (i) and (iv): Seven-index tensors, WL,wv 
^u/rrfa'rf' usec ^ t0 re P resent a gPEPS representation of the system's 
ground state wave functions for an infinite-size system, with s being 
a physical index, /, r, u, and d denoting the inner indices. Here, /', r', 
u' and d' are horizontal and vertical grading indices, respectively, (ii) 
and (v): The pictorial representation of a gPEPS \ifi) with horizontal 
and vertical grading bonds, which are used to absorb a two-site gate 
acting on horizontal and vertical bonds, respectively, (iii) and (vi): 
The unit cells of an infinite gPEPS with horizontal and vertical grad- 
ing bonds, respectively, made of four seven-index tensors W, X, Y, 
and Z. (vii) and (x): Double tensors wi u d r j', j and w tt ird; U 'd' we formed 
from the seven-index tensors W and theirs complex conjugates W* 
with horizontal and vertical grading bonds, respectively, (viii) and 
(xi): The tensor networks (TNs) for the norm of gPEPS's with hor- 
izontal and vertical grading bonds, respectively, (ix) and (xii): The 
unit cells of the TNs for the norm of gPEPS's with horizontal and 
vertical grading bonds, respectively. 



their complex conjugates (see Fig. [TJvii) and Fig. 02 X )X re- 
spectively. As such, one may visualize the norm for a gPEPS 
as a TN, as shown in Fig. [TJviii) and Fig.[Tfxi), with their unit 
cells plotted in Fig.QJix) and Fig. [TJxii). With the double ten- 
sors w, x, y and z as the building blocks, one may form the 
one-dimensional transfer matrix E\, which is a Matrix Prod- 
uct Operator on an infinite strip (see Fig. [3). The left and right 
eigenvectors corresponding to the largest eigenvalue of Ei 
are iMPS's, from which one may form the zero-dimensional 
transfer matrix Eq (see Fig. EJi)). The largest left and right 
eigenvectors of Eo are defined in Fig. SJii), which, together 
with those of E\, form the environment tensors. In addition, 
an auxiliary vector V R is defined by absorbing the tensors S3, 
S4, 1S 2 , El, y, and z, as visualized in Fig. EJiii). This enables 
us to compute the ground state energy for the XW bond, as 
shown in Fig.Uiv). 

The same procedure may be used to update the gPEPS ten- 
sors and to read out a physical observable for other bonds. 
However, different conventions should be adopted for hori- 
zontal and vertical bonds. 

We stress that the update procedure above is not optimal, 
in the sense that it does not produce the best approximate 




FIG. 2: (color online) The procedure to update the gPEPS tensors 
T w and T x and the singular value matrix A m . via absorbing the action 
of a two-site gate U„. (i) the two-site gate U a is applied onto the 
gPEPS. (ii) A single tensor © is formed by contracting the tensors 



rv, T x , A„, A z _ 



A yw , and the gate U a . (iii) Reshape the 



tensor © into a matrix M. (iv) A singular value decomposition (S VD) 
is performed for the matrix M, followed by a truncation, with only 
the D largest singular values retained in the updated singular matrix 
A n . w . (v) Reshape the matrices U and V into the tensors U and V. (vi) 
Recover the diagonal matrix A m , A zx , A xz , A ny , A yw , and update the 
tensors T w and T x to be T w and r„. 
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one - dimensional transfer matrix E, 

a. 




FIG. 3: (color online) The iMPS used to approximate the largest 
eigenvector of the one-dimensional transfer matrix E\ , shown here 
as an Matrix Product Operator on an infinite strip. Here, we need to 
absorb two-site nonunitary gate acting on an iMPS. The iMPS turns 
out to be the largest eigenvector of the transfer matrix E\ if A\, Ai, 
A-j, and A 4 converge after the transfer matrix E{ is acted on the iMPS 
enough times. 



gPEPS representation for each imaginary time slice during 
the imaginary time evolution. As such, our update proce- 
dure can only be used to produce the system's ground state 
wave functions, but not for real time evolution from a pre- 
scribed initial state. A similar situation occurs for an MPS 
algorithm |9|]. This drawback may be remedied if one uses 
the same strategy as the iPEPS algorithm 11711 . which is opti- 
mal in the above sense. That is, in order to absorb a two-site 
gate, one needs to compute the environment tensors, i.e., the 
left and right largest eigenvectors of both the one-dimensional 
and zero-dimensional transfer matrices for each time slices. 
Therefore, the update problem is reduced to a four-site sweep 
procedure that consists of successively solving a set of linear 
equations J3]. However, this requires to update the environ- 
ment tensors as we update the gPEPS tensors W, X, Y, and Z 
for each two-site gate, so it is much less efficient. 
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FIG. 4: (color online) The ground state energy per bond is com- 
puted by contracting the environment tensors, with the XW bond as 
an example, (i) The right and left largest eigenvectors of the trans- 
fer matrix E[ are denoted by tensors Ei, £2, S3, S4, and E 2 , 
E 4 , respectively. Here, the O-dimensional transfer matrix Eq is visu- 
alized, (ii) The largest left and right eigenvectors V L and V R of the 
one-dimensional transfer matrix E . (iii) An auxiliary vector V R is 
defined by absorbing the tensors X3, E4, 1.' 2 , ~L' 3 , y, and z- (iv) The 
ground state energy for the XW bond is computed by contracting a 
tensor T with the tensors V L , V R , Si, E2, Sj, and E^. Here, T is de- 
fined by the Hamiltonian density h X w an d the tensors x, w, and 
w*. 
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FIG. 5: (color online) For J > 0.95, there is a line of phase separation 
(PS). For J < 0.95, no PS occurs. Here, we have chosen D = 4. For 
J = 0.4 and n = 0.4968, denoted as A, the extended s-wave pairing 
order parameter < A >= 0.083 + 0.1 li, with the ground state energy 
per site e = -0.9498 for D = 6. For J = 3.0 and n = 0.1273, 
denoted as B, the extended s-wave pairing order parameter < A >= 
0.010 + 0.055i, with the ground state energy per site e = -0.4157 for 
D = 4. 



Simulation of the two-dimensional t — J model. We test the 
algorithm with the two-dimensional t - 7 model described by 
the Hamiltonian 111 80 : 



H = ~ t Yj mc ^ c i<r + H c ),p] + 3 Z (Si ' S J " \ ni (9) 



<lj>cr 



<'./> 



where S, are spin 1/2 operators at a lattice site i, V is the 
projection operator excluding double occupancy, and t and 7 
are, respectively, the hoping constant and anti-ferromagnetic 
coupling between the nearest neighbor sites < ij >. Hereafter, 
we shall choose t — 1 for brevity. 

At half filling (i.e., n = 1, with n being the number of elec- 
trons per site), the t - 7 model reduces to the two-dimensional 
Heisenberg model. In this case, the algorithm yields the 
ground state energy per site e - -1.16757, for the trunca- 
tion dimension D = 4, quite close to the QMC simulation 



result e = -1.16807 UJ, |20J. Away from half filling, the 
model exhibits different behaviors for small and large anti- 
ferromagnetic coupling 7, see Fig. For 7 > 0.95, there 
is a line of PS. For 7 < 0.95, no PS occurs. This agrees 
qualitatively with the results based on the high-temperature 
expansions [13]. Note that our result for the transition point 
J c = 3.43 at low electron density is quite close to the exact 
value J c - 3.4367 lElll . Here, we have chosen D = 4. 

In the homogeneous regime, it turns out that the algorithm 
does not yield much conclusive results, due to the fact that the 
truncation dimension D currently accessible is quite small (up 
to B = 6). However, signals of extended s-wave supercon- 
ductivity are observed in two regimes: the first is the regime 
for 2 < 7 < 3.43 at low electron density, and the second 
is the regime which starts at least from 7 = 0.4 at (almost) 
quarter filling. For 7 = 0.4 and n = 0.4968, denoted as 
A in Fig. [5] the extended s-wave pairing order parameter 

< A >= 0.083 + 0.11;, with the ground state energy per site 
e = -0.9498 for B = 6. For 7 = 3.0 and n = 0.1273, denoted 
as B in Fig. [5] the extended s-wave pairing order parameter 

< A >= 0.010 + 0.055/, with the ground state energy per site 
e - -0.4157 for B = 4. It remains unclear whether or not 
these two points are continuously connected. 

Given that the bottleneck of the algorithm to achieve higher 
precision is the smallness of the truncation dimension B, we 
expect that our data may be significantly improved for a larger 
truncation dimension B by incorporating the symmetry into 
the algorithm Il22ll . Indeed, even for the anti-ferromagnetic 
order parameter at half filling, the currently accessible trunca- 
tion dimensions are still too small. 

Summary and outlook. We have developed a numerical al- 
gorithm to find a gPEPS representation of the ground state 
wave functions for translationally invariant strongly correlated 
electronic systems on infinite-size lattices in two spatial di- 
mensions. It is tested for the two-dimensional t — J model 
at and away from half filling, with truncation dimensions up 
to 6. We are able to locate a line of PS, which qualitatively 
agrees with the results based on the high-temperature expan- 
sions [13J. It is proper to stress that the location of the line 
may vary if the truncation dimension is increased, although 
the variation might be small (especially at low electron den- 
sity), due to the fact that PS can be seen from a consideration 
based on energetics [ 23], whereas the ground state energy per 
site we computed is reasonable, compared to the exact values 
for a small (4 X 4) cluster. In addition, the model exhibits 
an extended s-wave superconductivity for 7 = 0.4? at quarter 
filling. However, we emphasize that the currently accessible 
truncation dimensions are not large enough, so it is necessary 
to incorporate the symmetry of the system into the algorithm, 
in order to achieve results with higher precision. This is cur- 
rently under investigation. 

After this work was completed, we have become aware 
of a preprint by T. Barthel, C. Pineda, and J. Eisert, 
arXiv:0907.3689, in which an alternative contraction scheme 
for the fermionic PEPS is discussed in the context of 
fermionic operator circuits. This work is supported in part 
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by the National Natural Science Foundation of China (Grant 
Nos: 10774197 and 10874252) and the Natural Science Foun- 
dation of Chongqing (Grant No: CSTC, 2008BC2023). 
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